Recurrent activity in neuronal avalanches

A new statistical analysis of large neuronal avalanches observed in mouse and rat brain tissues reveals a substantial degree of recurrent activity and cyclic patterns of activation not seen in smaller avalanches. To explain these observations, we adapted a model of structural weakening in materials. In this model, dynamical weakening of neuron firing thresholds closely replicates experimental avalanche size distributions, firing number distributions, and patterns of cyclic activity. This agreement between model and data suggests that a mechanism like dynamical weakening plays a key role in recurrent activity found in large neuronal avalanches. We expect these results to illuminate the causes and dynamics of large avalanches, like those seen in seizures.

are inter-avalanche firing times. Here we investigate that assumption. From the mouse data set we extracted the inter-fire times of each neuron and calculated both the average time between fires within the same avalanche and the average time between fires that were not in the same avalanche. In Figure 7 We plot the mean of each neuron below, with error bars showing the standard error of the mean. For each neuron, the average time between intra-avalanche fires (red) is orders of magnitudes smaller than the average time between inter-avalanche fires (green). This separation of the means was verified for each neuron with a hypothesis test that the two means were the drawn from the same normal distribution. We found the hypothesis to be rejected (p ∼ 0) for all neurons in the system. Thus, the assumption of these separated timescales has been justified for our data. We note that causal webs' do allow for simultaneous avalanches-yielding the possibility of a neuron having negative inter-avalanche firing times. We show in magenta the averages with these negatives values removed and see that the change to the results are negligible.

D: Synaptic time delay in simulations
Our model assumes that activity happening during an avalanche is much quicker than the time between consecutive avalanches i.e., that the neurons effectively communicate instantaneously. Here we use the model simulations to understand how our results change with the addition of a more biologically accurate synaptic time delay. In the original simulation protocol, the i th neuron firing at time t affected all other neurons j at time t+1. We show in Figure 8 the effect of delaying this until time t + τ for τ = 1, 2, 4, 20 simulation steps. Though we find that increasing time delay decreases the avalanche size cutoff, the effect of recurrent activations persists as an upward bending of the S vs. N curve.

E: Sub-sampling
The experimental setup suggests that roughly one quarter of the total number of neurons in the cortex slice are being recorded by the MEA. This is not because the electrode density is too sparse to see all of the neurons; in fact the electrode density is so high, all of the neurons above it can be identified with very high accuracy through triangulation [3]. The other three quarters of the neurons actually come from the other three quarters of the cortex slice that is not directly above the MEA. So, we have neurons outside of the view ('outside' neurons) and neurons inside the view ('recorded' neurons). Not recording the entire cortical slice allows the possibility of an 'outside' neuron firing and causing a 'recorded' neuron to fire as well.
To account for this in our model, we performed a simple sub-sampling procedure on our simulation output that allowed us to emulate this experimental detail; namely that the experiment involved an entire slice but only recorded from a quarter of it. If our experimental system consisted of N 'recorded' neurons, we would simply simulate 4N neurons, but only consider firing events from every fourth neuron in the system. This way, when we construct the avalanches, neurons are truly receiving information from 'outside' neurons, which causes them to fire seemingly unprovoked. We have thus emulated the experimental situation of neurons from outside the field of view influencing neurons within the field of view, some times referred to as a 'sub-sampling' effect.
The results in the main paper (Figure 3 a,b) are shown for simulations that had this sub-sampling procedure performed and we can see again consistency between the data and experiment. We have also verified that this sub-sampling was not affecting the presence of the recurrent activity. We show in Figure 9,a how the S vs. N curves change as we sub-sample. We find that the subsampling procedure enhances the recurrent activity in our simulations. We point out however that this enhancement through sub-sampling is not sufficient to account for the extent of the recurrent activity seen in the experiments alone ( Figure 9,b). We can reproduce experimental avalanche size CCDF's and S vs. N curves with our simulations with weakening alone or with weakening including sub-sampling, but not with sub-sampling alone. We also find the sub-sampling procedure can induce a clustered bi-test signature in simulations where it once was not present without sub-sampling, this trend is described in detail in the next section.

F: Bi-test for model simulations
Our simple slip model with weakening has been shown through large systemsize simulations to admit two unique signatures related to temporal correlations between avalanches [see chapter 2 in [5]]. Specifically, predictions about the intervals between consecutive avalanche start times can be made which admit a signature that differs between small and large avalanches. This signature, the bi-test of inter-avalanche start times (outlined in the manuscript), for small avalanches is either Poissonian when there is no weakening (Figure 10, Blue), or temporally clustered if there is weakening (Figure 10, Red). Though these signatures are seen in the bare MFT model, we can make the clustering vanish by making the width of the disorder distribution large enough (Figure 10, Green). Further, once we add inhibitory neurons, the clustering signature is almost certain to vanish (Figure 10, magenta). With both inhibitory neurons and high disorder width, as we need for our simulation to match the model, we also fail to recover clustered small avalanche inter-start times (Figure 10, cyan).
Since the clustered signature seen in the data is not predicted by our simulations (i.e., to be a result of weakening), what causes it? We have investigated the effects of sub-sampling the simulations to more closely replicate the experiments (see Section E above). We find that, under sub-sampling of our simulations, the clustered bi-test signature for small avalanches is recovered (Figure 11). We were not able to get a sub-sampled simulation to all at once match the size distribution, S vs. N and the Bi-test, but we propose that sub-sampling is a possible reason for the appearance of the clustering signature in the data. We leave a perfect implementation of a sub-sampling scheme to future work, for this paper, we simply show that it is a possible explanation.

G: Disorder Configurations
Each time we run the simulation, a new set of arrest stresses is drawn from a uniform distribution. This means that each simulation will not have exactly the same evolution over time. Specifically, the stress configuration of the system can stay relatively uniform for the duration of the simulation, or it can evolve into a bunched state; see [2]. This can result in a handful of large events (significantly larger than the rest) obscuring the fit to experiment (Figure 12, disorder realization # ′ s 1,3-6,8-10). Importantly the statistical features we are concerned with, for instance, the power-law slope of the avalanche size CCDF, do not change with these different realizations of disorder ( Figure 12). With that being said, the experiments have one realization of disorder, which we are apriori unaware of. For this reason, in the main paper ( Figure. 4), we highlight one disorder configuration that compares best with the experiment (disorder realization #2, Figure 12).  Bending in S vs. n. Fire signals of cultured cortex recorded with a 512-electrode array and processed into causal webs were used to extract from each avalanche the number of unique neurons (n) and size (S). The two are plotted vs. each other above on a log-log scale, with the black line S(n)=n, shown here for five critical mouse data sets (a-e) and one critical rat data set (f) not shown in the main paper. Fire signals of cultured cortex recorded with a 512-electrode array and processed into causal webs were used to extract the distribution of the mean number of times each neuron fires in an avalanche, for avalanches within specified size. They are shown here for five critical mouse data sets (a-e) and one critical rat data set (f) not shown in the main paper. Fire signals of cultured cortex recorded with a 512-electrode array and processed into causal webs were used to extract cycles of recurrent activity shown here for five critical mouse data sets (a-e) and one critical rat data set (f) not shown in the main paper.

Fig. S 5:
Bi-test. Fire signals of cultured cortex recorded with a 512-electrode array and processed into causal webs were used to extract avalanche intereventtimes. The cumulative distribution function of H i = ∆t i /(∆t i + ∆τ i /2) is shown here for large (empty markers) and small (filled markers) avalanches. ∆t i is taken as the minimal interevent time, min(|t i − t i−1 |, |t i+1 − t i |), and τ i is the subsequent interevent time (|t i−1 −t i−2 | or |t i+2 −t i+1 |). Poisson distributed random interevent times fall directly on the dashed line, while perfectly periodic interevent times fall on the dotted line. The green regime is a signature of temporal clustering while in the purple regime is a signature of temporal quasiperiodicity. They are shown here for five critical mouse data sets (a-e, blue, square) and one critical rat data set (f, red, triangle) not shown in the main paper.

Fig. S 6:
Traditional avalanche definitions over-emphasize recurrent firing. Fire signals of cultured mouse cortex recorded with a 512-electrode array processed two ways: (blue) sorting fires into time bins of width δt, defining an avalanche as any consecutive set of non-empty bins bounded by empty bins [1]. (black) sorting fires into causally connected space (neuron)-time pairs, defining an avalanche as any string of causally connected fires [6]. Avalanche sizes are plotted versus the number of unique neurons that fired in the avalanche.   Simulations of the mouse system (N=434) were ran with W=1.9, ϵ e = 0.5 ϵ i = 0.7, γ = 4.8. (Left) simulations with nonzero weakening show the subsampled system has good agreement with experiment (red line). (Right) the sub-sampled system without weakening does not replicate the data nearly as well as the simulations with weakening (red line). From this figure we conclude that sub-sampling alone cannot account for the recurrent activity seen in the data. Simulations of the mouse system (N=434) were run with low disorder (W=0.1), ϵ e = 0.5 ϵ i = 0.7, γ = 4.8. With out any weakening, the small avalanches show a near Poissonian bi-test signature (blue). With the addition of weakening, the signature resembles clustering (red). However, the addition of strong disorder (W=1.9, green) or inhibitory neurons (magenta) or both (cyan) leads to the destruction of this clustered signature. From this investigation we provide a possible explanation for why the clustered signature is not present in our model simulations that include inhibitory neurons and large disorder. Simulations of the mouse system (N=434) were run with W=1.9, ϵ e = 0.5 ϵ i = 0.7, γ = 4.8. The solid line shows the full simulation (with 4N neurons) and the dashed line that same simulation after being sub-sampled by a factor of 4. From this figure we conclude that sub-sampling could be partially responsible for the clustered signature we see for small avalanches in the data. We show the CCDF of avalanche size for 10 different realizations of disorder (or 'trials') in varying colors. The experimental trial from the main paper, mouse, is shown in purple and the disorder averaged CCDF is shown in black. The left figure shows the CCDF with out sub-sampling, the right figure shows the CCDF with sub-sampling. All configurations on the left match the tail of the distribution fairly well, along with the disorder average. There is discrepancy for smaller sized avalanches though, and this discrepancy seems to relax with sub-sampling (right figure). The trade off is that there are sometimes large avalanches in the sub-sampled (and accordingly quadrupled in size) system, as explained in the SI, Sec G. Still, an individual simulation trial (trial 2 and 7, for instance) compare quite well to one experimental trial.
where ∆t i is taken as the minimal interevent time, min(|t i − t i−1 |, |t i+1 − t i |), and τ i is the subsequent interevent time (|t i−1 − t i−2 | or |t i+2 − t i+1 |) shown here for example Poissonian (left), periodic (center) and quasi-periodic (right) fire-trains. For applications of the Bi-test to plasticity see [4]. , where N t is the total number of neurons and a factor of -1 is added when neuron i and j are not of the same type. Avalanche sizes are calculated, along with total number of participating neurons, shown here plotted against each other on a loglog scale. Fire signals of cultured cortex recorded with a 512-electrode array and processed into causal webs were used to extract from each avalanche the number of unique neurons (n) and size (S) shown in black circles. (Top) Non-symmetric coupling: The excitatory to inhibitory weights were reduced by a factor of one (red dots) one half (blue dots) and one third (green dots). (Bottom) Weakening needed for bending in model: The model simulations are plotted with (red dots) and without (blue dots) dynamical weakening.